% fig3.m（部分修正）
clear; clc;

% 无量纲化参数
param.C1 = 1; param.L1 = 1; param.R1 = 0.8; param.E1 = 0.7;
param.C2 = 1; param.L2 = 1; param.R2 = 0.8; param.E2 = 0.7;
param.LA = 10; param.RS = 10; 
param.a = 1.0; param.V0 = 1.0;

% 设置洛伦兹系统参数并存入结构体
param.sigma = 10; 
param.rho = 28; 
param.beta = 8/3;

% 初始条件（u0, v0, w0）
X0 = [20; 0; 0]; 

% 时间区间
tspan = [0, 1000];

% RK4 求解洛伦兹系统
h = 0.1;  % 步长
[T, X] = rk4(@(t, X) lorenz(t, X, param), tspan, X0, h);

% 提取 u, v, w
u = X(1, :);
v = X(2, :);
w = X(3, :);

% 外部声信号
A = 0.9; f = 0.05; 
param.vPCFunc = @(t) A * cos(2*pi*f*t);

% 初始条件
X0 = [0.2; 0.1; 0.2; 0.1; 0]; 

% 调用RK4求解器
[T, Y] = rk4(@(t, X) inductive(t, X, param), tspan, X0, h);

% (a) 子图 (a)
% subplot(5,1,1);
M = 0;  % 设置 M=0
param.iChaosFunc = @(t) M * interp1(linspace(0, 1000, length(u)), u, t);
[T, Y] = rk4(@(t, X) inductive(t, X, param), tspan, X0, h);
plot(T, Y(1,:)/param.V0, 'k'); 
hold on;
plot(T, param.vPCFunc(T)/param.V0, 'r');
ylim([-2, 2]);
title('(a) M=0');

% % (b) 子图 (b)
% subplot(5,1,2);
% M = 0.03;  % 设置 M=0.03
% param.iChaosFunc = @(t) M * interp1(linspace(0, 1000, length(u)), u, t);
% [T, Y] = rk4(@(t, X) inductive(t, X, param), tspan, X0, h);
% plot(T, Y(1,:)/param.V0, 'k'); 
% hold on;
% plot(T, param.vPCFunc(T)/param.V0, 'r');
% ylim([-2, 2]);
% title('(b) M=0.03');
% 
% % (c) 子图 (c)
% subplot(5,1,3);
% M = 0.05;  % 设置 M=0.05
% param.iChaosFunc = @(t) M * interp1(linspace(0, 1000, length(u)), u, t);
% [T, Y] = rk4(@(t, X) inductive(t, X, param), tspan, X0, h);
% plot(T, Y(1,:)/param.V0, 'k'); 
% hold on;
% plot(T, param.vPCFunc(T)/param.V0, 'r');
% ylim([-2, 2]);
% title('(c) M=0.05');
% 
% % (d) 子图 (d)
% subplot(5,1,4);
% M = 0.1;  % 设置 M=0.1
% param.iChaosFunc = @(t) M * interp1(linspace(0, 1000, length(u)), u, t);
% [T, Y] = rk4(@(t, X) inductive(t, X, param), tspan, X0, h);
% plot(T, Y(1,:)/param.V0, 'k'); 
% hold on;
% plot(T, param.vPCFunc(T)/param.V0, 'r');
% ylim([-2, 2]);
% title('(d) M=0.1');
% 
% % (e) 子图 (e)
% subplot(5,1,5);
% M = 0.3;  % 设置 M=0.3
% param.iChaosFunc = @(t) M * interp1(linspace(0, 1000, length(u)), u, t);
% [T, Y] = rk4(@(t, X) inductive(t, X, param), tspan, X0, h);
% plot(T, Y(1,:)/param.V0, 'k'); 
% hold on;
% plot(T, param.vPCFunc(T)/param.V0, 'r');
% ylim([-2, 2]);
% title('(e) M=0.3');
